In [1]:
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
import seaborn as sns
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
import warnings
warnings.simplefilter('ignore')
from sklearn.datasets import load_boston
In [2]:
boston_dataset = load_boston()
In [3]:
type(boston_dataset)
Out[3]:
In [4]:
boston_dataset
Out[4]:
In [5]:
dir(boston_dataset)
Out[5]:
In [6]:
print(boston_dataset.DESCR)
In [7]:
print(boston_dataset.data)
In [8]:
print(boston_dataset.feature_names)
In [9]:
print(boston_dataset.filename)
In [10]:
print(boston_dataset.target)
In [11]:
data=pd.DataFrame(data = boston_dataset.data, columns = boston_dataset.feature_names)
data
Out[11]:
In [12]:
data['PRICE'] = boston_dataset.target
In [13]:
data
Out[13]:
In [14]:
data.head()
Out[14]:
In [15]:
data.tail()
Out[15]:
In [16]:
data.count()
Out[16]:
In [17]:
pd.isnull(data)
Out[17]:
In [18]:
pd.isnull(data).any()
Out[18]:
In [19]:
data.info()
In [20]:
# plotting...visualisation
In [21]:
plt.hist(data['PRICE'])
plt.show()
In [22]:
plt.hist(data['PRICE'], bins = 30)
plt.show()
In [23]:
plt.figure(figsize =(18, 8))
plt.hist(data['PRICE'], bins = 30, ec = 'black')
plt.xlabel('Price in 1000\'s of Dollars', fontsize = 20)
plt.ylabel('Nr of Houses', fontsize = 20)
plt.show()
In [24]:
data['PRICE'].mean()
Out[24]:
In [25]:
data['PRICE'].median()
Out[25]:
In [26]:
# mean and median are different.....SO that means the data is not normally distributed
In [27]:
# what is a normal distribution
In [28]:
# seaborn............
In [29]:
sns.distplot(data['PRICE'])
plt.show()
In [30]:
plt.figure(figsize =(18, 8))
sns.distplot(data['PRICE'], hist=True,kde=False, hist_kws={'ec':'black', 'color': 'red', 'alpha': 0.3})
plt.show()
In [31]:
plt.figure(figsize =(18, 8))
sns.distplot(data['PRICE'], hist=True,kde=True, hist_kws={'ec':'black', 'color': 'red', 'alpha': 0.3})
plt.show()
In [32]:
plt.figure(figsize =(18, 8))
sns.distplot(data['PRICE'], hist=True,kde=False, hist_kws={'ec':'black', 'color': 'red', 'alpha': 0.3})
plt.show()
In [33]:
plt.figure(figsize =(18, 8))
sns.distplot(data['PRICE'], hist=False,kde=True, kde_kws={'color': 'red'})
plt.show()
In [34]:
data['RM'].mean()
Out[34]:
In [35]:
data['RM'].median()
Out[35]:
In [36]:
plt.figure(figsize =(18, 8))
sns.distplot(data['RM'], hist=True,kde=False, hist_kws={'ec':'brown', 'color': 'black', 'alpha': 0.3})
plt.show()
In [37]:
data['CHAS']
Out[37]:
In [38]:
data['CHAS'].value_counts()
Out[38]:
In [39]:
data['CHAS'].values
Out[39]:
In [40]:
data['CHAS'].index
Out[40]:
In [41]:
type(data['CHAS'])
Out[41]:
In [42]:
plt.figure(figsize =(10,6))
plt.hist(data['RAD'], ec = 'black', color = 'navy', alpha = 0.5)
plt.xlabel('Accessibility to Highways')
plt.ylabel('Nr of Houses')
plt.show()
In [43]:
accessibility=data['RAD'].value_counts()
accessibility
Out[43]:
In [44]:
type(accessibility)
Out[44]:
In [45]:
accessibility.index
Out[45]:
In [46]:
type(accessibility.index)
Out[46]:
In [47]:
accessibility.values
Out[47]:
In [48]:
type(accessibility.values)
Out[48]:
In [49]:
plt.figure(figsize =(10,6))
plt.bar(accessibility.index, accessibility.values, ec = 'black', width = 0.5)
plt.show()
In [50]:
data.min()
Out[50]:
In [51]:
data.max()
Out[51]:
In [52]:
data.median()
Out[52]:
In [53]:
data.mean()
Out[53]:
In [54]:
data.describe()
Out[54]:
In [55]:
# why plotting.......
Anscombe's quartet comprises four data sets that have nearly identical simple descriptive statistics, yet have very different distributions and appear very different when graphed. Each dataset consists of eleven (x,y) points. They were constructed in 1973 by the statistician Francis Anscombe to demonstrate both the importance of graphing data before analyzing it and the effect of outliers and other influential observations on statistical properties.
In [56]:
# Correlation
In [57]:
data['PRICE'].corr(data['RM'])
Out[57]:
In [58]:
plt.figure(figsize = (18, 12))
plt.scatter(data['PRICE'], data['RM'], color = 'red', alpha = 0.4)
plt.show()
In [59]:
data['PRICE'].corr(data['PTRATIO'])
Out[59]:
In [60]:
# why is this negative?
In [61]:
# I will give you 10 minutes for the answer....
In [62]:
plt.figure(figsize = (18, 12))
plt.scatter( data['PTRATIO'],data['PRICE'], color = 'navy', alpha = 0.4)
plt.show()
In [63]:
data['NOX'].corr(data['PRICE'])
Out[63]:
In [64]:
plt.figure(figsize = (18, 12))
plt.scatter( data['NOX'],data['PRICE'], color = 'black', alpha = 0.4)
plt.show()
In [65]:
data['CHAS'].corr(data['PRICE']) # it means the houses close to the river banks are priced marginally higher
Out[65]:
In [66]:
data.corr()
Out[66]:
In [67]:
# Heat map
In [68]:
plt.figure(figsize =(16,10))
sns.heatmap(data.corr(), annot = True, annot_kws = {'size':14}, fmt = '.1g', mask=None)
plt.show()
In [69]:
# mask : boolean array or DataFrame, optional
# If passed, data will not be shown in cells where ``mask`` is True.
# Cells with missing values are automatically masked.
In [70]:
data.corr().shape
Out[70]:
In [71]:
masker = np.zeros_like(data.corr())
In [72]:
masker.shape
Out[72]:
In [73]:
masker
Out[73]:
In [74]:
triangle_indices=np.triu_indices_from(masker)
masker[triangle_indices] = True
masker
Out[74]:
In [75]:
plt.figure(figsize =(16,10))
sns.heatmap(data.corr(), annot = True, annot_kws = {'size':14}, fmt = '.1g', mask=masker, linewidths = 0.3)
plt.xticks(fontsize = 10)
plt.yticks(fontsize = 10)
plt.show()
In [76]:
# Home work....
In [77]:
# why is there a -ve correlation between the PRICE and LSTAT
# why is there a +ve correlation between the TAX and NOX
# why is there a -ve correlation between the DIS and NOX
In [78]:
# why is there a -ve correlation between the PRICE and LSTAT
# LSTAT % lower status of the population
# Lower status....low income...not very educated....more crimes...more fights...more problems in the area...
# would you want to go and buy a house there?.....no...not at all...
# hence houses in such areas will have low demand..and hence the negative correlation...between PRICE and LSTAT
# to prove that there is a corelation between...LSTAT and CRIM lets examine the data...
# heat map....for the CORR
# there is a positive correlation between CRIM and LSTAT...that proves the point..
In [79]:
# why is there a +ve correlation between the TAX and NOX
# NOX.....nitrous oxide pollutant found in areas where there are industries...
# TAX.....looks to be higher in areas where there are more industries...that sounds logical..
# INDUS proportion of non-retail business acres per town?....what is that?.....is that Industrial Area or Residential Area?
# if INDUS is industrial Area...then there must a positive correlation between INDUS and NOX....INDUS and TAX
# we have found that there is a positive correlation between INDUS and TAX...0.7
# we have found that there is a positive correlation between INDUS and NOX....0.8
# so INDUSTRIAL AREAS HAVE TO PAY HIGHER TAX, ALSO THEY ARE MORE POLLUTED..
In [80]:
# why is there a -ve correlation between the DIS and NOX...
# DIS weighted distances to five Boston employment centres...
# NOX nitrous oxide concentration..
# Employment centres were mostly factories those days...located in industrial areas..
# far away from the industrial area then there will less NOX...and hence the neagative correalation...
# there must be a -ve correlation between DIS and TAX..?
# can we check? ...its negative -0.5....
In [81]:
# seaborn plotting...
In [82]:
# import seaborn as sns
In [83]:
# kind : { "scatter" | "reg" | "resid" | "kde" | "hex" }, optional
In [84]:
sns.set()
sns.set_style('darkgrid')
sns.jointplot(x =data['NOX'] , y =data['DIS'], color='red', alpha = 0.3, height=15, ratio=6, kind='scatter')
plt.show()
In [85]:
sns.set()
sns.set_style('white')
sns.lmplot(x = 'TAX', y = 'RAD', data = data, height=10, aspect=2)
plt.show()
In [86]:
sns.set()
sns.set_style('white')
sns.lmplot(x = 'RM', y = 'PRICE', data = data, height=10, aspect=2)
plt.show()
In [87]:
sns.set()
sns.set_style('whitegrid')
sns.lmplot(x = 'RM', y = 'PRICE', data = data, height=10, aspect=2, ci = None)
plt.show()
In [88]:
# Multivariable regression...
In [89]:
data
Out[89]:
In [90]:
prices = data['PRICE']
In [91]:
prices
Out[91]:
In [92]:
features = data.drop('PRICE', axis =1)
In [93]:
features
Out[93]:
In [94]:
# split this data into train and test...
# train...x_train and y_train...x_train will hold 80% of the features...y_train will hold 80% the prices
# test....x_test and y_test....x_test will hold 20% the featues...y_test will hold 20% the prices
# from sklearn.model_selection import train_test_split
In [95]:
X_train, X_test, y_train, y_test = train_test_split(features, prices, test_size = 0.2, random_state = 10)
In [96]:
X_train
Out[96]:
In [97]:
y_train
Out[97]:
In [98]:
X_test
Out[98]:
In [99]:
y_test
Out[99]:
In [100]:
len(features)
Out[100]:
In [101]:
len(X_train)
Out[101]:
In [102]:
round(len(X_train)/len(features),2)*100 # train_size
Out[102]:
In [103]:
round(len(X_test)/len(features),2)*100 # test_size
Out[103]:
In [104]:
regr = LinearRegression() # training the model
regr.fit(X_train, y_train)
Out[104]:
In [105]:
regr.coef_
Out[105]:
In [106]:
regr.intercept_
Out[106]:
In [107]:
pd.DataFrame(data =regr.coef_ , index = X_train.columns, columns = ['Coef'])
Out[107]:
In [108]:
regr.score(X_train, y_train)
Out[108]:
In [109]:
regr.score(X_test, y_test)
Out[109]:
In [110]:
data['PRICE']
Out[110]:
In [111]:
plt.figure(figsize= (16,8))
plt.hist(data['PRICE'], bins = 50)
plt.show()
In [112]:
data['PRICE'].skew()
Out[112]:
In [113]:
# Data Transformation...
In [114]:
A = 90
B = 50
In [115]:
A-B
Out[115]:
In [116]:
np.log(A)
Out[116]:
In [117]:
np.log(B)
Out[117]:
In [118]:
np.log(A)-np.log(B)
Out[118]:
In [119]:
y_log=np.log(data['PRICE'])
y_log
Out[119]:
In [120]:
y_log.skew()
Out[120]:
In [121]:
plt.figure(figsize= (16,8))
plt.hist(y_log, bins = 50)
plt.show()
In [122]:
plt.figure(figsize= (16,8))
plt.subplot(1,2,1)
plt.hist(data['PRICE'], bins = 50, color = 'red', alpha = 0.4)
plt.subplot(1,2,2)
plt.hist(y_log, bins = 50, alpha = 0.5)
plt.show()
In [123]:
# do a regression after data transformation ...that means...prices will be log prices..
In [124]:
data
Out[124]:
In [125]:
log_prices = np.log(data['PRICE'])
In [126]:
features = data.drop(['PRICE'], axis = 1)
In [127]:
log_prices
Out[127]:
In [128]:
features
Out[128]:
In [129]:
X_train, X_test, y_train, y_test =train_test_split(features, log_prices, test_size = 0.2, random_state = 10)
In [130]:
regr.fit(X_train, y_train)
Out[130]:
In [131]:
regr.coef_
Out[131]:
In [132]:
pd.DataFrame(data=regr.coef_, index=X_train.columns, columns=['coef']) # after data transformation...
Out[132]:
In [133]:
regr.intercept_
Out[133]:
In [134]:
regr.score(X_train, y_train)
Out[134]:
In [135]:
regr.score(X_test, y_test)
Out[135]:
In [136]:
# So, data transformation has improved the scores....
# So we will keep the transformation..
In [137]:
# Plots with sns....
In [138]:
sns.set()
plt.figure(figsize=(20,6))
sns.distplot(y_log)
plt.title(f'Log price: The skew is {round(y_log.skew(),2)}', fontsize = 25)# note the f string used to print the y_log.skew()
plt.xlabel('LOG PRICES ', fontsize =14, color = 'green' )
plt.show()
In [139]:
skew=round(data['PRICE'].skew(),2)
sns.set()
plt.figure(figsize=(20,6))
sns.distplot(data['PRICE'])
plt.title(f'Actual price: The skew is {skew}', fontsize = 25)
plt.xlabel( 'ACTUAL PRICES ', fontsize =14, color = 'green' )
plt.show()
In [140]:
# we are going to conduct further tests on the data.....
# do we need all these features ...`13 of them..?
# Is there a problem of multicollinearity amongst the features...check VIF....?
# Are all these features really related to the price..target....check for p values???
In [141]:
# hypothesis testing...
In [142]:
# EXP HYP:-...... Varma is in love with Rani.....
# NULL HYP:-....Varma is not in love with Rani....
# GET DATA:-.... based on the null hypothesis...Varma is not in love with Rani.....
# then we are calculating the Prob...of what?... Varma is not in love with Rani....
# if P(Varma is not in love with Rani) > 0.05....null hpy may be possible...and varma might not be in love with Rani..
# if p(Varma is not in love with Rani) < 0.05 ...then null is not possible..So maybe Varma is in love with Rani...
In [143]:
# sklearn.....
# statsmodels....will allow me to do a lot of tests..also do the regression..
In [144]:
results=sm.OLS(y_train, X_train).fit()
# we have not added the constant to the exog...so we did not get the constant..intercept and also got wrong values of the coef..
results.params
Out[144]:
In [145]:
X_train
Out[145]:
In [146]:
sm.add_constant(X_train)
Out[146]:
In [147]:
results=sm.OLS(y_train, sm.add_constant(X_train)).fit() # used statsmodels...SM
results.params
Out[147]:
In [150]:
round(results.pvalues,3)
Out[150]:
In [151]:
pd.DataFrame({'coef': results.params, 'p-value': round(results.pvalues, 3)})
Out[151]:
In [152]:
# Multicollinearity....
# VIF....if VIF > 10...there is multicollinearity...then we might have to drop the feature...just might have to...
$$ TAX = \alpha _0 + \alpha _1 RM + \alpha _2 NOX + ... + \alpha _{12}LSTAT $$¶
$$ VIF _{TAX} = \frac{1}{(1 - R _{TAX} ^ 2)} $$¶
In [154]:
variance_inflation_factor(exog =sm.add_constant(X_train).values, exog_idx = 10 )
Out[154]:
In [155]:
variance_inflation_factor(exog =sm.add_constant(X_train).values, exog_idx = 1)
Out[155]:
In [160]:
for i in range(len(sm.add_constant(X_train).columns)):
print(variance_inflation_factor(exog =sm.add_constant(X_train).values, exog_idx = i),'\n')
In [161]:
# INDUS and AGE have p values greater than 0.05
# VIF looks okay....
In [163]:
# Bayesian information criterion.....BIC
# BIC?....Deals with what?....Deals with overfitiing or overlearning...
In [164]:
results.bic
Out[164]:
In [165]:
# Build our models.....
# All features and the log prices....Model 1
In [166]:
# without INDUS and the log prices....Model2
In [167]:
# Without INDUS and AGE with log prices...model3
In [168]:
# we will look at the BIC of all these models...
# we will look at the r sqaure of all these models..
In [170]:
X_incl_const = sm.add_constant(X_train)
model = sm.OLS(y_train, X_incl_const)
results_1 = model.fit()
In [171]:
results_1.params
Out[171]:
In [172]:
results_1.bic
Out[172]:
In [173]:
results_1.rsquared
Out[173]:
In [174]:
X_incl_const
Out[174]:
In [175]:
X_incl_const = X_incl_const.drop(['INDUS'], axis =1)
model = sm.OLS(y_train, X_incl_const)
results_2 = model.fit()
In [176]:
results_2.bic
Out[176]:
In [178]:
results_2.rsquared
Out[178]:
In [180]:
X_incl_const = X_incl_const.drop(['AGE'], axis =1) # already removed INDUS
model = sm.OLS(y_train, X_incl_const)
results_3 = model.fit()
In [181]:
results_3.bic
Out[181]:
In [182]:
results_3.rsquared
Out[182]:
In [183]:
# So what?....the model to go with is model...3...with the least BIC...
In [ ]: